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Abstract — Cell-centered finite-volume (CCFV) schemes have certain attractive properties for 
the solution of the equations governing compressible fluid flow. Among others, they provide a 
natural vehicle for specifying flux conditions at the boundaries of the physical domain. Unfor- 
tunately, they lead to slow convergence for numerical programs utilizing them. In this report a 
method for investigating and improving the convergence of CCFV schemes is presented, which 
focuses on the effect of the numerical boundary conditions. The key to the method is the compu- 
tation of the spectral radius of the iteration matrix of the entire discretized system of equations, 
not just of the interior point scheme or the boundary conditions. 

1 Introduction 

The traditional finite- volume schemes are based on cell-centered schemes. In these schemes the flow 
quantities are stored at the center of the computational cell and the metrics defining the shape 
of the cell are at the cell vertices or the nodes of the computational grid. An alternative method 
is the cell- vertex based scheme, where the physical and the geometric quantities are stored at the 
cell vertices. In both methods the viscous flow equations are solved using the integral form of 
the Navier-Stokes equations, which guarantees that the conservation principles of the governing 
equations are also satisfied in the discretized form. 

Cell-centered schemes have been in wide-spread use for a number of years [5, 6, 9], but more 
recently cell-vertex schemes have been gaining popularity [4, 8, 10, 11]. Cell-centered finite-volume 
schemes have certain desirable properties for the solution of compressible fluid flow problems. They 
are generally more accurate than cell-vertex based schemes, use fewer operations per point since 
averaging is done over smaller stencils, and also lead to more physical boundary conditions, since flux 
computations at the boundary of the computational domain coincide with actual grid boundaries. In 
contrast, cell-vertex schemes usually employ extrapolations from the interior of the grid to produce 
boundary fluxes. In addition, cell-vertex schemes require an extra step to recover the update of the 
physical quantities at the cell vertices. 

The accuracy claim made above for the cell-centered scheme is disputed by Roe [10] and Rossow 
[11], who maintain that for non-smooth meshes that are typical of realistic geometries cell-vertex 
schemes are at least first order accurate, whereas cell-centered schemes are of zeroth order. In other 



words, the truncation error remains finite for vanishing grid spacing. However, it has been shown 
[13, 14] that even for non-smooth meshes, cell-centered schemes are always first order if the fluxes 
at the cell faces are computed properly. On smooth meshes cell-centered schemes are second order 
accurate, but cell-vertex schemes are second order only under special circumstances. 

One reason why cell-centered schemes are still not much used is that they tend to lead to 
very slow convergence in viscous flow situations with explicit boundary value updates. This has 
not been a major difficulty in the earlier application of these schemes when the flow solvers were 
explicit. However, with the advent of implicit methods such as the alternating-direction implicit 
schemes [1] developed for the Navier-Stokes equations, the slow convergence rates have been a major 
handicap. The slow convergence rates of cell-centered schemes were first noted when some of the 
popular finite-difference codes, such as CNS [2] and ARC2D [12], were converted to cell-centered 
finite- volume codes (see for example the papers by Klopfer [7] and Yoon [15]). The time steps 
had to be reduced by more than an order of magnitude to maintain numerical stability. A similar 
deterioration of convergence was noted with cell-centered unstructured grid schemes [3]. As will be 
shown in this study, the major cause of the slow convergence is the explicit application of certain 
boundary conditions. Finite-difference and cell-vertex finite-volume schemes are impervious to this 
problem. 

The slow convergence can be relieved by imposing the boundary conditions implicitly. However 
implicit boundary condition implementations are often difficult or impossible to achieve with popular 
structured-grid solution algorithms such as Beam- Warming ADI schemes [1] or LUSGS [15], and 
with unstructured-grid algorithms [3] as well. Here we will demonstrate, through the study of 
a succession of model equations of increasing complexity, that the judicious choice of relaxation 
factors for the explicit update of boundary values can render such updates as efficient numerically 
as would be a completely implicit implementation. 


2 Steady, scalar 1-D model equation 

A simple model equation for viscous fluid flow is the linearized, steady, one-dimensional convec- 
tion/diffusion equation. 


du 




d 2 u 


dx dx 2 ^ 

Assuming unit density throughout, the convection speed c and viscosity (i can be combined into the 
Reynolds number: Re = f c//z. We also introduce the cell Reynolds number on a uniform grid with 

mesh spacing A x : Re c = f Re A x . Applying central differencing to equation 1 leads to the following 
discrete equation. 


(1 -I- Re c /2)iij_i — 2iq + (1 — Re c /2)u i+ i = 0 . (2) 

A second-order-accurate finite-volume formulation of the Dirichlet boundary conditions u\ x=a = u a , 
u\ X -b = Ub is obtained by setting 

\{u\ +u 2 ) = u a , i(u n _i +u n ) = u b , ( 3 ) 
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where the true physical boundaries at x = a and x — b are at half-way points x a and Xf, between 
interior grid points and ghost points. See Figure 1. 


u a “b 



Figure 1: One-dimensional finite-volume grid 


If the one-dimensional problem is embedded in an operator-split multi-dimensional problem, 
boundary conditions schemes are often applied separately from the interior-point advancement. 
Such schemes are usually implemented by first updating the ghost-point values U\ and u n using 
Equation 3 to ensure that the boundary conditions are satisfied, meanwhile keeping the interior- 
point values constant. Subsequently, the interior-point values are updated (in this case by using 
Equation 2), but now the ghost-point values are kept unchanged. If we use the tilde symbol ‘"’to 
indicate the intermediate ghost-point values after application of the boundary conditions, we can 
write the above scheme as a two-step semi-implicit procedure: 

Step 1: 



Combining Equations 4 and 5, we obtain a single-step iterative scheme whose iteration matrix N 
is given by: 
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The eigenvalues of N determine the convergence rate of the scheme. These are found by setting 
Nq = Aq: 


-92 = A<7i , 

0 = A{(1 + Re c /2)qi-i - 2$ + (1 - Re c /2)q i+1 } , * = 2, ...,n - 1 , (7) 

9n— 1 = Xq n • 

The solution to this set consists of two distinct families. The first one has n — 2 eigenvectors 
that all correspond to an eigenvalue of zero. This eigenspace is spanned by the n — 2 basis vectors: 
span{ei, e 3 , e 4 , . . . , e n _ 3 , e„_ 2 , e n }. The second — non-trivial — family has nonzero eigenvalues, which 
determine the actual convergence of the scheme. The general solution to the equations for the 
eigenvectors at interior grid points is of the form 

J i Re 

Qi = a + Pr ' , with r = ~ L ' ( 8 ) 

1 2 

The constants a and /3 are determined by applying the two ‘boundary conditions’ of Equation 7, 
which are written in matrix form: 


{ A + 1 Ar + t 2 
\X + 1 At" + r"- 1 



(9) 


This equation only has nontrivial solutions if the determinant of the matrix on the left hand side is 
zero, which yields the following eigenvalues: 


A - I r - rn ~ 2 


Weak convection 

In case of weak convection we can write r ~ 1 + Re c and r k « 1 + kRe c , so 

A 1>2 = -l,^ + C>(Jfe*). (11) 

n — 1 

Apparently, the effect of convection on the rate of convergence of the scheme vanishes (to first order) 
for small values of the cell Reynolds number. Since n > 3 for a finite-volume scheme, A 2 is always 
negative, and smaller in absolute value than 1. But A 2 approaches —1 for large n, and Ai = — 1, so 
the method will never converge, although the interior points are updated using an implicit scheme. 
The lack of convergence is due to the boundary-value update scheme. The fact that the eigenvalues 
of the iterative scheme are both (close to) —1 suggests, however, that underrelaxation may help. 
Consequently, we replace the boundary value update Equation 4 by: 



4 



All other steps in the determination of the eigenvalues of the iteration matrix stay the same, and 
we now find: 


Ai, 2 « 


( 20 - 1 ), 


(20- l)n + 3 


Minimizing the spectral radius of the iteration matrix yields: 


3 — 2n 


p(N) 


3 — 2n 


Obviously, excellent convergence is obtained with this choice of the underrelaxation factor 9. For 
large grid sizes the value of 1/2 for 9 is completely satisfactory. 


Strong convection 

The situation changes rather dramatically in case the convection is strong (large Re c ). We introduce 
the small parameter e = Ij Re c . Again, the first non-trivial eigenvalue of A is —1. But the value of 
A 2 depends on the parity of n: 

A 2 = 1 — 8 (n — 2)e 2 -f 0(e 3 ) for n even (15) 

= — — - + 0(e 2 ) for n odd (16) 

n — 1 

Regardless of the parity of n, however, A 2 approaches 1 for large grid sizes as the cell Reynolds 
number goes to infinity, which indicates that underrelaxation will not speed up convergence. Indeed, 
if we underrelax the boundary conditions, we find: 

Ai ,2 w (20 — 1), lor ^ — y(l — 0) + 0. (17) 

It follows that A 2 is always close to or equal to 1, independent of the value of 9. 

The situation may be remedied by using a more physically relevant boundary condition at the 
downstream boundary of the domain. Assuming a positive convection speed, we replace the Dirichlet 
boundary condition at x = b by a Neumann condition. A second-order-accurate implementation of 
the condition %z\ x =b = 0 is obtained by setting: 


= <- 


(18) 


Underrelaxing the Neumann and Dirichlet conditions using 9 n and Op, respectively, produces the 
following system: 
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We repeat the eigenvalue/eigenvector analysis carried out above. The components of the non- 
trivial eigenvector q can again be written as = a + /?r\ but now the definition of r is changed 
to: r = ( 2e — l)/(2e + 1). Notice that now lim e | 0 T = -1 (large cell Reynolds number), whereas 
before it was 1 (small cell Reynolds number). The small (2 x 2)-determinant equation producing 
the non-trivial eigenvalues of the iteration matrix (cf. Equation 9) becomes: 

A-l A T-T 2 + 9 n {t 2 -t) / 2 qn 

A + 1 — 2#o At” + r n_1 — 0 D {r n + r" -1 ) ' v ’ 

The expressions for the roots Ai, 2 are not very insightful in general. We investigate the case where 
the cell Reynolds number approaches infinity, so r approaches —1. 

For an even number of grid points (n), Equation 20 reduces to: 


2A 2 — (Op + $at)A + 1 + (1 — 29d){1 — 26 n) — 0 . 


( 21 ) 


The roots of this equation are both zero (and hence produce a spectral radius equal to zero) for 
0 D = ±l/\/2 , 0 N = —0 D . Notice, however, that a negative value of a relaxation factor corresponds 
to an overrelaxation (i.e. an extrapolation) of a boundary value, which is prone to instability. If we 
limit the relaxation factors to positive values and set 6 n equal to zero, we find an optimal spectral 
radius of p(N) = \/2 — l(s=s 0.4) by choosing 9p — 2(\/5 — 1). 

Unfortunately, the picture is not so rosy in case the number of grid points is odd. Equation 20 
becomes singularly perturbed for large cell Reynolds numbers (coefficient of A 2 disappears), and we 
have to take limits carefully. Whereas the root corresponding to the reduced equation (Re c = oo) 
is given by 

, l-(l-20*)(l-2 0 D ) , ooX 

A - 4-2(0 w + «„) ’ ' ' 

which can always be made to equal zero, the other root is asymptotically proportional to the cell 
Reynolds number, irrespective of the choice for the relaxation factors; the scheme always diverges 
for large Re c . 

The solution to the problem — not too surprisingly — is to use upwind differencing in case of 
strong convection. The resulting iteration matrix for the Neumann/Dirichlet condition is: 


( 1 

I 1 + Rg c 


N = 


V 


—2 — Re c 
1 Re c 


\ ( 9d 6d — 1 

0 


—2 - Re r 1 


1 ) 


\ 


° 

1 ~ 0N &N ) 


(23) 


Equation 20 is still valid, provided r is defined by 1 + Re c . For large cell Reynolds numbers it 
reduces to: 


(20d — A — l)(0;v — A) = 0 , 


(24) 
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so a spectral radius of approximately zero can be obtained by setting Qd = 1/2 and On — 0. For 
small cell Reynolds numbers the eigenvalue equation is: 

(n - 1)A 2 + (20 D (1 -n) + 2- 0 N ) A + (2 0 D - l)(n - 2 + 0 N ) - 1 = 0 . (25) 

For large n (many grid points) the spectral radius has the form: p{N) = 1 — the convergence 

deteriorates for vanishing convection and large grid sizes, and the Neumann outflow boundary 
condition should not be used. Reverting to Dirichlet conditions, underrelaxed by 50%, restores the 
excellent convergence also obtained when using central differencing. 

We now combine the previous analyses into a rule of thumb: 

Apply Dirichlet boundary conditions on inflow boundaries and always underrelax by 50%. In case 
of strong convection, use upwind differencing and apply Neumann boundary conditions on outflow 
boundaries, but do not underrelax these conditions. 


3 Time-dependent, scalar 1-D model equation 

We now investigate the convergence properties of finite- volume formulations of the linearized, time- 
dependent, one-dimensional convection/diffusion equation: 


du 

dt 


du 


d 2 u 


+ c di~^dz* 


(26) 


We apply the Euler-implicit method for the time integration (^-superscript) and central differencing 
in space (z-subscript) at all interior grid points. 

(1 + Re c l2)u k ffl - (2 + v)u k+l + (1 - Re c /2)u k +{ = -uu k , with (27) 


def cAf 
CFL = — — 

Ax 


def 
V — 


Re r 


Ax 2 


CFL pAt 


Again, boundary values are updated as in Equation 4, where k now signifies the time step instead 
of the iteration number. Advancement of the interior-point values takes place through solution of 
Equation 27. The resulting two-step scheme can again be combined into a single equation, whose 


N — 


matrix is 

given by 

(cf. 

Equation 6): 



i 



\ 

-1 

0 -1 

\ 

i + f 

—2 — v 

1 - 

Re 

2 


—v 



1 + Ss. 
~ 2 

-2 

- V 1 - f 



—v 




1 J 


\ 

-1 0/ 


(28) 


Explicit expressions for the eigenvalues of N cannot be obtained in general, but some limiting 
cases are easily identified. 


lim JV = 1 , lirnN = 1 . 

v— H-oo 14 0 


(29) 
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Whereas the first limit is expected and acceptable — for infinitely small time steps convergence will 
be infinitely slow as well — the second is highly undesirable, since it implies that very large time steps 
will lead to very slow convergence as well. This is an artifact of the application of the boundary 
conditions of the finite-volume scheme, and we attempt to solve the convergence problem in the 
same vein as before. 

Weak convection 

In Figure 2 the effect of underrelaxing Dirichlet boundary conditions for a weakly convective case 
(Re c — 2) on a small grid (n = 40) is depicted. The cell Reynolds number is chosen such that 
the discretization matrix for the interior point scheme is guaranteed to be diagonally dominant for 
all finite sizes of the time step, a property often judged desirable in fluid flow computations. The 
dashed and solid lines indicate the amplification factor in the absence of relaxation and with 50% 
underrelaxation, respectively. For curiosity’s sake we also plot the amplification factor corresponding 
to 30% underrelaxation (dotted line). 
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Figure 2: Amplification factor for underre- Figure 3: Amplification factor for under- 

laxed Dirichlet boundary conditions; weak relaxed Dirichlet boundary conditions; no 

convection (Re c = 2) convection (Re c — 0) 

Notice that all three curves collapse onto the same branch for larger values of v, i.e. as the 
time step goes to zero. The improvement through the 50% underrelaxation for small values of v 
is apparent. The fact that the amplification factor does not decrease strictly monotonically as v 
approaches zero is due to the presence of convection. As the cell Reynolds number vanishes, the 
underrelaxed amplification factor becomes strictly monotone and reduces to zero for infinite time 
step (see Figure 3). The curve corresponding to the nonrelaxed boundary conditions, however, 
remains non-monotone. 

If we underrelax the mildly convective case between 50% and 30% , the intersection point with 
the vertical axis (amplification factor for infinite time step) does not change much, although the 
minimum amplification factor for nonzero v may decrease. For relaxation factors outside the above 
range the intersection point actually moves up, and gains are diminished. 

It is also apparent from Figures 2 and 3 that even underrelaxation by 50% does not completely 
save the difference scheme from deterioration in the case of total absence of convection. Only 
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for infinite time step does the amplification factor go to zero, but for finite time step the factor 
approaches 1 ever more rapidly as the cell Reynolds number drops. 

Strong convection 

As in the steady analysis, we apply upwind differencing and Neumann boundary conditions at the 
outflow boundary for the strongly convective case, and modify Equation 28 accordingly. Now it 
is most useful to study the behavior of the difference scheme for a certain fixed CFL number. We 
first investigate the case of infinitely strong convection. The outflow boundary condition cannot be 
enforced anymore since the spatial differencing reduces in the limit to a one-sided difference, so the 
set of equations determining the eigenvalues of the amplification matrix becomes (cf. Equation 7): 

-x 2 = \X\ , 

— Xj/CFL = A (*,_! - (1 + l/CFL)x.) . (30) 

Solving the recursive interior-point equation and substituting the boundary condition yields the 
following asymptotic result for large CFL numbers: p(N ) = 1 — 2/CFL. Applying underrelaxation 
of the Dirichlet (upstream) boundary condition gives for large values of the CFL number: 



Re c tec 

Figure 4: Amplification factor for underre- Figure 5: Amplification factor for un- 

laxed Dirichlet boundary conditions; small derrelaxed Dirichlet boundary conditions; 

CFL number (CFL = 1) medium CFL number (CFL = 10) 

Obviously, 0 D = 1/2 is again a good choice. For cell Reynolds numbers of finite magnitude 
we apply the same underrelaxation of the Dirichlet boundary condition, and leave the Neumann 
boundary condition again untouched. The results for varying values of the CFL number are show in 
Figures 4, 5, and 6, together with the amplification factors for no underrelaxation of the Dirichlet 
condition (dashed line), and underrelaxation of 30%. Although 0d = 1/2, On = 0 is not always the 
optimal choice, it is always good, asymptotically best (large CFL), and easy to apply. 
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Figure 6: Amplification factor for underrelaxed Dirichlet boundary conditions; large CFL number (CFL = 
1000) 


4 Unsteady, scalar 2-D model equation 


We now extend the analysis to a two-dimensional, linear model equation for the convection/diffusion 
process. 


du du du ( d 2 u d 2 u\ 

dt Cx dx + ^ dy ^ \<9a: 2 + dy 2 ) 


(32) 


This problem is solved on a square mesh of (n — 2) x (n — 2) interior points The grid is depicted in 
Figure 7, including the ghost points (crosses) that lie outside the grid boundaries. We also employ 
the ghost points on the corners of the exterior grid (open squares) to arrive at a square total number 
of points. But since they never enter the discretization scheme or the solution, we keep the values at 
these corner ghost points fixed. We apply central differencing to all spatial derivatives (5 Xi denotes 
the central difference operator in the X{ direction), and backward differencing (Euler Implicit) to 
the time derivative to arrive at the following discretization. 


I + CFL X 




+ CFLj, 



u k+l — 14* +1 . 


(33) 


Here we have defined: Re Cx . = and CFL Xi = c ^ t . If one or both of the convection speeds 

c Xi are zero, equation 33 is simplified accordingly. As in the one-dimensional finite-volume case, we 
update ghost point values such that the average of every couple of values at grid points straddling 
the boundary equals the Dirichlet boundary value, or symbolically: u,+u x = Ub, where the subscript 
b stands for ‘boundary’. Application of the ghost point value update operator £2, which is the direct 
analogue of the one-dimensional operator, gives the following time advancement scheme. 


(I + L x + Ly)u k+1 = Qu k . 


(34) 


Note the introduction of L x and L y as shorthand for the spatial discretization operators. From 
equation 34 we derive two additional schemes through ADI approximate factorization. Under the 
assumption that L x and L y are small compared to I, a reasonable approximation is: 

(/ + L x ) (I + Ly)u k+l = Qu k . (35) 
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Figure 7: Two-dimensional finite-volume grid 


For large At, however, this assumption is not valid, and we have to adopt the so-called delta 
form of the ADI method. This is derived from equation 33 by introducing the update quantity 
A u k — u k+l — f ~lu k , and by applying the approximate factorization only to the implicit part (‘left 
hand side’) of the operator that multiplies A u k . The resulting delta form can be combined with 
the non-delta form into a single equation. 

(I + L X )(I + L y )u k+l = (/ + (!- a)L x L y )Vtu k . (36) 

If a equals zero we obtain the delta form, whose final solution is independent of the size of the time 
step. For any other a , the solution at convergence depends on the size of the time step, due to the 
factorization error, which is proportional to aL x L y . Algebraic analysis of the spectral radius of the 
iteration matrix corresponding to equation 36 is too complex, so we resort to numerical investigation 
instead. Figure 8 show's the spectral radius of the iteration matrix on a grid of 40 x 40 points for the 
non-delta form of the ADI algorithm, and 9 for the delta form. Since the direct-inversion scheme 
does not incur a factorization error, its delta and non-delta forms are identical. From Figure 8 
we deduce that the non-delta forms benefit significantly from underrelaxation of the boundary 
conditions for large time steps (small v). In fact, the ADI scheme converges as rapidly as the 
direct-inversion scheme. But if we insist on a solution independent of the size of the time step, the 
delta form of the ADI algorithm must be used, w’hich gains virtually nothing from underrelaxation 
(see Figure 9, where the curves for ADI with and without boundary update relaxation coincide). 
Since direct inversion is not practical for many calculations of engineering value, we must find ways 
of speeding up the ADI algorithm without giving up solution accuracy at convergence. 

One approach to accomplishing this is to let the switch parameter a in the hybrid scheme 
represented by equation 36 slide from one to zero during the course of the iteration process. We use 
the fact that the convergence factor for the underrelaxed hybrid scheme behaves approximately as 
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Figure 8: Amplification factor for non- 
delta form (a = 1); no convection 
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Figure 9: Amplification factor for delta 
form (a = 0); no convection 


follows for large time steps and for small a. 

p(a) «l-sa, (37) 

where s is a positive constant. This relationship is depicted in Figure 10, where the computed 
spectral radius for the iteration matrix belonging to the 50% underrelaxed ADI scheme on the 
40 x 40 grid is shown as a function of the switch parameter. The non-dimensional time step is set at 
an arbitrarily large value of 1000. Note that keeping a constant at some value other than zero will 
always speed up the convergence, but at the cost of final solution accuracy. We now analyze the 
convergence of the entire iterative process for some choice of the way in which a goes to zero, i.e. 
lim^^oo a k = 0, where k is the number of the iteration or time step. Since the convergence is faster 



Figure 10: Amplification factor for hybrid ADI scheme; no convection 


for values of a away from 0, we select a fairly slowly decreasing function, namely a — zj (fc+1), where 
z is some positive constant. Without loss of generality, we can set z equal to the computationally 
convenient 1/s. Consequently, we obtain p k d = p(a k ) « 1 - l/{k + 1). The ‘convergence’ factor for 
p steps of the hybrid iterative scheme is: 


n pk 



i 

p + 1 


(38) 
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Evidently, the scheme is converging very slowly, since after p steps the compound convergence factor 
is only l/(p + 1). The error of the scheme, defined by e p = u°° — u p , is estimated as follows. At 
convergence, the solution u°° satisfies equation 34, i.e., 

{I + L x + L y )u°° = Qu°°. (39) 


Subtracting equation 36 from equation 39 yields the following expression: 


e 


p+i 


[(/ + L X )(I + Ly)}-' [(/ + (1 - a p )L x L y ) Qe p + 
a p [(I + L X )(I + Ly)]- 1 L x L y nu°° 

R(a p )e p + a p t°° 

n + ( 52 

1 \/c = l 2 = 1 / 


(40) 


where R is the iteration matrix of the scheme defined by equation 36, and t°° is a quantity dependent 
only on the converged solution. The norm of the error can be bounded as follows: 


ie p+i n < n p* ii e ‘n + (evhi n a) ii*° 

k = \ \ A: — 1 i= 1 J 


(41) 


Consequently, the solution error generally decreases at least as slowly as indicated by the compound 
convergence factor (equation 38). For the choice of a = 1 /{s(k + 1)), we find 


r +1 u < 


He'll f 1 l|t“ll 

p+ 1 k 2 s 


which leads to the limiting result: 


lim ||e p+1 || < 

p -> OO " M 



(42) 


(43) 


This means that the scheme never converges to the correct solution, which appears to contradict 
the fact that the compound convergence factor does approach zero.. But the computed convergence 
factors apply to an ever-changing discretization scheme, due to the ‘hybridization 1 , and no conclu- 
sions can be drawn about the eventual vanishing of their product per se. We note that even in the 
case of a fairly slowly changing a, the scheme stalls. If we make a approach 0 more slowly, the final 
solution error will be even larger. If more quickly, then convergence will be even slower. It thus 
follows that the hybrid scheme does not have a fundamentally better convergence behavior than 
the plain delta form of the ADI scheme, and relaxation of the boundary condition application will 
not improve its performance. 
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5 Conclusions 


We have shown that it is possible to overcome the problems associated with the application of time- 
lagged boundary value updates for cell-centered finite-volume schemes for simple one-dimensional 
scalar, linear model equations. The appropriate technique relies on underrelaxing the boundary 
value updates to annihilate or reduce oscillations. In two dimensions the same technique works 
well for an approximately factored scheme, but the factorization error renders the solution useless 
for large time steps. When the delta form of the scheme is used to reduce the factorization error, 
convergence becomes very slow, even for a diffusion-dominated flow, and using underrelaxation of 
boundary value updates is ineffective. 
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